#include "gauss-legendre.h"

/* Valores iniciais para o método de Gauss-Legendre */
void gl_init_values(mpf_t a, mpf_t b, mpf_t t, mpf_t p)
{
    /* Calcula a0 = 1 */
    mpf_init_set_ui(a, 1);

    /* Calcula b0 = 1/sqrt(2) */
    mpf_init_set_ui(b, 2); mpf_sqrt(b, b); mpf_ui_div(b, 1, b);

    /* Calcula t0 = 1/4 */
    mpf_init_set_d(t, 0.25);

    /* Calcula p0 = 1 */
    mpf_init_set_ui(p, 1);
}

/* Inicializa estrutura paralela */
void gl_init_param(gl_param *p)
{
    /* Inicia com os valores iniciais */
    gl_init_values(p->a[0], p->b[0], p->t[0], p->p[0]);

    mpf_inits(p->a[1], p->a[2], p->b[1], p->b[2], p->t[1], p->p[1], p->p[2], p->pi[0], p->pi[1], NULL);

    p->stop = B_FALSE;
}

/* Limpa estrutura paralela */
void gl_clear_param(gl_param *p)
{
    mpf_clears(p->a[0], p->a[1], p->a[2], p->b[0], p->b[1], p->b[2], p->t[0], p->t[1],
            p->p[0], p->p[1], p->p[2], p->pi[0], p->pi[1], NULL);
}

/* Realiza o cálculo de a */
void gl_a_value(mpf_t a1, mpf_t a0, mpf_t b0)
{
    mpf_add(a1, a0, b0); mpf_div_ui(a1, a1, 2); // a1 = (a0+b0)/2
}

/* Realiza o cálculo de b */
void gl_b_value(mpf_t b1, mpf_t a0, mpf_t b0)
{
    mpf_mul(b1, a0, b0); mpf_sqrt(b1, b1); // b1 = sqrt(a0*b0)
}

/* Realiza o cálculo de t */
void gl_t_value(mpf_t t1, mpf_t t0, mpf_t p0, mpf_t a0, mpf_t a1)
{
    mpf_sub(t1, a0, a1); mpf_mul(t1, t1, t1); mpf_mul(t1, t1, p0); mpf_sub(t1, t0, t1); // t1 = t0-p0*(a0-a1)^2
}

/* Realiza o cálculo de p */
void gl_p_value(mpf_t p1, mpf_t p0)
{
    mpf_mul_ui(p1, p0, 2); // p1 = 2*p0
}

/* Realiza o cálculo de Pi */
void gl_pi_value(mpf_t pi, mpf_t a, mpf_t b, mpf_t t)
{
    mpf_add(pi, a, b); mpf_mul(pi, pi, pi); mpf_div(pi, pi, t); mpf_div_ui(pi, pi, 4); // pi = (a1+b1)^2/(4*t1)
}

/* Método de Gauss-Legendre sequêncial */
void gauss_legendre(mpf_t pi, long *it)
{
    mpf_t a[2], b[2], t[2], p[2];
    mpf_t old_pi;
    int i, i0, i1;

    /* Inicia com os valores iniciais */
    gl_init_values(a[0], b[0], t[0], p[0]);

    mpf_inits(a[1], b[1], t[1], p[1], old_pi, NULL);

    for(i = 0; ; i++)
    {
        i0 = INDEX2(i, 0);
        i1 = INDEX2(i, 1);

        /* Calcula os parâmetros */
        gl_a_value(a[i1], a[i0], b[i0]);
        gl_b_value(b[i1], a[i0], b[i0]);
        gl_t_value(t[i1], t[i0], p[i0], a[i0], a[i1]);
        gl_p_value(p[i1], p[i0]);

        mpf_set(old_pi, pi);
        
        /* Calcula o valor de pi */
        gl_pi_value(pi, a[i1], b[i1], t[i1]);

        /* Verifica convergência */
        if(has_converged(old_pi, pi)) break;
    }

    mpf_clears(a[0], a[1], b[0], b[1], t[0], t[1], p[0], p[1], old_pi, NULL);

    *it = i;
}

/* Método de Gauss-Legendre concorrente */
void gauss_legendre_conc(mpf_t pi, long *it)
{
    pthread_t thread[5];
    gl_param p;
    int i, j;

    gl_init_param(&p);

    for(i = 0; !p.stop; i++)
    {
        p.i = i;

        /* Cria uma thread para cada parâmetro */
        pthread_create(&thread[0], NULL, glt_calculate_a, &p);
        pthread_create(&thread[1], NULL, glt_calculate_b, &p);
        pthread_create(&thread[2], NULL, glt_calculate_t, &p);
        pthread_create(&thread[3], NULL, glt_calculate_p, &p);
        pthread_create(&thread[4], NULL, glt_calculate_pi, &p);

        /* Espera todas as threads acabar para sincronizá-las */
        for(j = 0; j < 5; j++)
            pthread_join(thread[j], NULL);
    }

    *it = i; /* Número de iterações */
    mpf_set(pi, p.pi[INDEX2(i, 1)]);
    gl_clear_param(&p);
}

/* Gauss-Legendre Thread para cálculo de a */
void *glt_calculate_a(void *arg)
{
    gl_param *p = (gl_param *) arg;
    int i0 = INDEX3(p->i, 0);
    int i1 = INDEX3(p->i, 1);

    gl_a_value(p->a[i1], p->a[i0], p->b[i0]);

    return NULL;
}

/* Gauss-Legendre Thread para cálculo de b */
void *glt_calculate_b(void *arg)
{
    gl_param *p = (gl_param *) arg;
    int i0 = INDEX3(p->i, 0);
    int i1 = INDEX3(p->i, 1);

    gl_b_value(p->b[i1], p->a[i0], p->b[i0]);

    return NULL;
}

/* Gauss-Legendre Thread para cálculo de t */
void *glt_calculate_t(void *arg)
{
    gl_param *p = (gl_param *) arg;
    int i0, i1, j0, j1;

    /* "Bolhas" para cálculo de t */
    if(p->i < 1) return NULL;

    i0 = INDEX3(p->i-1, 0);
    i1 = INDEX3(p->i-1, 1);
    j0 = INDEX2(p->i-1, 0);
    j1 = INDEX2(p->i-1, 1);

    gl_t_value(p->t[j1], p->t[j0], p->p[i0], p->a[i0], p->a[i1]);

    return NULL;
}

/* Gauss-Legendre Thread para cálculo de p */
void *glt_calculate_p(void *arg)
{
    gl_param *p = (gl_param *) arg;
    int i0 = INDEX3(p->i, 0);
    int i1 = INDEX3(p->i, 1);

    gl_p_value(p->p[i1], p->p[i0]);

    return NULL;
}

/* Gauss-Legendre Thread para cálculo de Pi */
void *glt_calculate_pi(void *arg)
{
    gl_param *p = (gl_param *) arg;
    int j0, i1, j1;

    /* "Bolhas" para cálculo do Pi */
    if(p->i < 2) return NULL;

    i1 = INDEX3(p->i-2, 1);
    j0 = INDEX2(p->i-2, 0);
    j1 = INDEX2(p->i-2, 1);

    gl_pi_value(p->pi[j1], p->a[i1], p->b[i1], p->t[j1]);

    /* Se houve convergência, a variável stop é setada como true */
    if(has_converged(p->pi[j0], p->pi[j1])) p->stop = B_TRUE;

    return NULL;
}

